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High-precision evaluation of Wigner’s d-matrix by exact diagonalization 
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The precise calculations of the Wigner’s d-matrix are important in various research fields. Due to the presence 
of large numbers, direct calculations of the matrix using the Wigner’s formula suffer from loss of precision. We 
present a simple method to avoid this problem by expanding the d-matrix into a complex Fourier series and 
calculate the Fourier coefficients by exactly diagonalizing the angular-momentum operator Jy in the eigenbasis 
of J-. This method allows us to compute the d-matrix and its various derivatives for spins up to a few thousand. 
The precision of the d-matrix from our method is about 10“'^' for spins up to 100. 


PACS numbers: 02.60.-x, 03.65.Fd, 02.20.-a, 21.60.-n 

I. INTRODUCTION 

The spin is a fundamental quantum object and an impor¬ 
tant candidate for various quantum technologies such as mag¬ 
netic resonance spectroscopy, quantum metrology, and quan¬ 
tum information processing. An essential requirement in 
these developments is the precise control of many spins or 
alternatively a large spin composed of the constituent spins. 
The simplest case of such control is the rotation around a 
fixed axis. Accurately describing this process requires high- 
precision calculations of the Wigner’s d-matrix M that 
quantifies the rotation of angle 9 around the y axis: = 

{j, m\ exp(-i9Jy)\j, n) 6 R, where |y, m) is an eigenstate of 
with eigenvalue m, i.e., Jf\j, m) - m\j, m). Hereinafter, H = I 
and i = V^. 

High-precision calculations of the d-matrix is of interest 
in quantum metrology iHl. For instance, let us consider 
an atomic Ramsey (or equivalently, Mach-Zehnder) interfer¬ 
ometer fed with all spins down as the paradigmatic setup of 
interferometric phase estimation. These spins then undergo 
an unknown phase shift 6 via the evolution exp{-i9Jy) inside 
the interferometer. Finally, by detecting the population imbal¬ 
ance at the output port of the interferometer, i.e., the mea¬ 
surement with respect to the output state exp(-i9Jy)\j, -j), 
one can record (2] + 1) possible outcomes. The outcome m 
occurs with probability Pm{6) = \{j, m\ exp{-i6Jy)\j, -j)'^ - 
[d'^ conditioned on the unknown parameter 6, thus 6 

can be inferred from appropriate data processing on the out¬ 
comes. This process, however, requires accurate evaluation 
of Wigner’s d-matrix. In addition, the ultimate sensitiviw of 
this estimation is determined by the Fisher information O-Q]: 
F(0) = Yjtn[dPm(ff)186]^IPm(6), which requires accurate eval¬ 
uation of the first-order derivative of Wigner’s d-matrix. 

In addition to quantum metrology, the Wigner’s d-matrix is 
closely related to spherical harmonics and Le gend re polyno¬ 
mials and is of interest in many other fields lISl- ll^ . However, 
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the calculation of the d-matrix for large spins (y » 1) suffers 
from a serious loss of precision, due to the presence of large 
numbers that exceed the floating-point precision in Wigner’s 
original formula M. To avoid this problem, the d-matrix 
has been calculated by means of recurrence relations ^ . This 
method still encounters severe numerical instability in the 
case of high spin, although a few remedies have been pro¬ 
posed IfSl - l^ . Recently, Gumerov and Duraiswami ll24l] have 
developed a new recursion relation for each subspace of spins, 
which greatly improves the stability. However, the achievable 
precision (i.e., the maximum absolute error) of their results re¬ 
mains unclear. Most recently, Tajima ll2^ proposed Fourier- 
series expansion of the Wigner’s d-matrix and convert the ac¬ 
curate evaluation of the d-matrix to that of the Fourier coef¬ 
ficients. Such a Fourier-series representation has been shown 
to be more useful in improving the numerical stability and 
the precision. However, each Fourier coefficient is still a sum 
of many large numbers that exceed the floating-point preci¬ 
sion, so it has to be evaluated with the assistance of a formula- 
manipulation software ll^ . 

In this paper, we put forward a very simple method to 
resolve the above large-number problem in evaluating the 
Fourier coefficients of Wigner’s d-matrix ll2^ . The essen¬ 
tial idea is to express these coefficients via the inner products 
{j, m\j,ff)y, where the eigenstates of Jy, i.e., {\j,ji)y} constitute 
an orthonormalized and completed set. To evaluate such inner 
products, we write down 7,. as a Hermitian matrix in the eigen¬ 
basis of J^. Then we numerically diagonalize the 7,. matrix to 
obtain the eigenstates and the inner products. Due to the nor¬ 
malization of \j,iu}y, the norm of each Fourier coefficient is 
not larger than unity, thus we avoid the large-number prob¬ 
lem in floating-point calculations. This method allows us to 
evaluate accurately the d-matrix and its various derivatives for 
much larger spins up to a few thousands, with an absolute er¬ 
ror (9(10 for the d-matrix and (9(/l0 for its kth-order 
derivative. 


II. FOURIER SERIES OF WIGNER’S D-MATRIX 

For large spin j, numerical calculation of the Wigner’s for¬ 
mula is subject to intolerable numerical errors because it is a 
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sum of many large numbers with alternating signs ffl: 



0 -. 2 j- 2 k+n~m, 0 -. 2 k+m-n 


, V 2) 


( 1 ) 


where 


where = e ' 2 '^'e are eigen¬ 

states of Jy and they constitute an ortho-normalized and com¬ 
pleted set, i.e., = 5^,^, and Y,f,\j,n)yy{j,iJ-\ = 1. 

Hereafter, we use \j, m) for the eigenstates of and |y,p)j, for 
the eigenstates of Jy. Comparing Eq. (|2|l and Eq. (IHi, we iden¬ 
tify the Eourier coefficients in Eq. (|2]l as 


Umn) (-1)^^'”-" VO' + m)\{j - m)\{j + n)\{j - n)\ 

^ (j - m - k)\(j + n - k)\(k + m - n)[k[ 

and k e [max(0, n - m), min(y - m, j + n)]. Taking Q(7r/2) 
as an example, the term k = y/2 has a very large magni¬ 
tude a 2V i, exceeding the floating-point precision 

when j » 1. 

To avoid this problem, Tajima 12^ has proposed to expand 
the Wigner’s d-matrix into a complex Eourier series: 


din,n{ff} = ^ 


-ifiG .{j,m,n) 
a 


( 2 ) 


This representation of the d-matrix is very useful and free 
from the large-number errors since each Eourier coefficient is 
less than or equal to 1 (see below). However, an accurate eval¬ 
uation of the Eourier coefficients by means of Eq. O 

remains nontrivial. This is because the large-number problem 
still exists in the coefficients: 

rain(;-mj'+n) 

^ wf"''''hi,{2j,2k + m-n), (3) 

k=ma\{0,n-m) 


where we have introduced an integral 

1 p27T , 0y2j-~X , ny\ 

V2J.A).-J^ (cos-) (sin-) C'-Sie 


>^0 

minUj+/y} 


i E 


(- 1 ) 


/-A/2 


/=max{0,-j+/i+A} 


2j - A \(A' 


j + fi-ijyij’ 


with (j) = A\l[l\{A - /)!]. When y » 1, some terms in the 
integral are still hug e (e.g., the term A - j, I - 0, and yu = 
-jl2). Tajima 112^ bypassed this problem by employing a 
symbolic computation software and then reducing the results 
to double-precision floating numbers. 


III. METHOD OF EXACT DIAGONALIZATION 

Here, instead of using Eq. (I3]l, we present a very simple 
method to calculate the Eourier coefficients that free 

from the above mentioned large-number problem. The key 
observation is that the d-matrix can be rewritten as 

dL,ni0) = {j, \j, n) = 2] Mj, lAyyij, ^^\i, n), (4) 

/^--j 


^U,m,n) ^ n)yy{j, ^i\i, n) = e'5(" (^) (^) ’ (5) 

which obeys the sum rule = {j, m\j, n) - 

Erom Eq.®, one can note that all the Eourier coefficients and 
hence the d-matrix for arbitrary 0 depend on di,^{0 - njl) 
only IHIlilll]. The elements of the d-matrix are real and 
obey the following symmetry properties (see Refs. fllH, and 
also Fig.®: 


dU{0) = di,n{-0) = = {-\r-’^diU0\ 

dUn- 0 ) = (-1)^'^'"<_„(0) = {-\y-"dL„,^ni 0 \ 
dL,n(7T + fl) = (-l)^'-"<_„(0) = (-iy^’"di^^„(0X 


from which one can easily obtain 

and as observed recently 

by Tajima j^ . 



FIG. 1: (Color online) Computed results of the Wigner’s d-matrix 
dm,ii(0) (with total spin j = 40) against m and n for 6 = n/6 (a), jijA 
(b), and ;r/2 (c), respectively. The dashed lines in bottom panel are 
the boundary of the central region, determined by Eq. ®. Outside 
the region, the values of are almost vanishing. 


Most importantly, the first result of Eq. © provides a very 
simple but accurate method to calculate the Fourier coeffi¬ 
cients and hence the Wigner’s d-matrix by solving the in¬ 
ner product {j,m\j,fi)y. To this end, we first express Jy as a 
(2 j -I- l)-dimensional Hermitian matrix: 

0 -X^jyl ] 

Xj 0 


0 -Xj 

X^j^i 0 j 
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where the matrix elements are determined by 
{j,m\Jy\j,n) ^ with the 

term X„, - yj( j + m)(j - m + 1) satisfying X±m - Xq:„,+i and 
X-j - 0. Next, we diagonalize the Hermitian matrix using 
standard numerical methods, e.g., the ZHBEV subroutine of 
LAPACK, or the EVCHF package of the IMSL, to obtain 
all the eigenvectors {\j,fi)y} and their probability amplitudes 
{j, m\j,fj.)y. Eor the simplest case j - 1/2, the matrix is indeed 
the transpose of the Pauli matrix o-y over two (i.e., cr^ 12), 
which gives two eigenvectors |l/2,+l/2)j,, corresponding 
to the eigenvalues +1/2. The inner products {j,m\j,jj.)y for 
j = 1/2 are given by <1/2,1/2|1/2,±1/2X. = 1/V2 and 
<l/2,-l/2|l/2,±l/2)y = +// V2. Inserting them intoEq. (01, 
one can obtain the Wigner’s d-matrix. 

The exact-diagonalization method has two advantages. 
First, the magnitude of {j,m\j,fj.)y and hence all the coeffi¬ 
cients in Eq. (|5]l are not larger than unity, since all the 

eigenvectors {\j,fJ.}y} are normalized. This provides a solu¬ 
tion to the large-number problem in Eqs. ([TJ and (l3]l. Second, 
the matrix in Eq. (|6ll is tridiagonal and Hermitian, which can 
be easily diagonalized. By diagonalizing the Jy matrix only 
once, all the Fourier coefficients and hence all the elements 
of the d-matrix can be obtained with given j and 6 (see the 
Supplemental Material ll^ i. This method allows us to cal¬ 
culate Wigner’s d-matrix for j up to a few thousands. More 
important, one can compute arbitrary k-th order derivative of 
the d-matrix with little additional cost, due to 

(7) 

This is indeed a common advantage of the Fourier-series rep¬ 
resentation. As comparison, one can find that the direct eval¬ 
uation of the first result of Eq. (|3 costs double even for the 
case k - \ since {j, m\ txp{—i6Jy){—iJy)\j, n) gives 

ddl,„(6) 1 r / /I 

, ( 8 ) 

which depends on two elements of the d-matrix. When n - 
± j, only one of them remains. 


IV. NUMERICAL RESULTS AND ERROR ANALYSIS 


Given the exact value t/gx and the numerically calculated 
value t/gomp of a matrix element d{„ „(6), the absolute error is 
defined as Aabs = l^fcomp - d^^l and the relative error is defined 
as Ai-gi = Kt/comp - <7ex)/<7exl- The exact values of the d-matrix 
elements are obtained by MATHEMATICA 10.0. Figure |2] 
shows that the absolute error Aabs ~ 10“'"^ even for a relatively 
larger spin j = 100 (see also Fig. |3ll. The lower panel of Fig. 
shows the relative error < 10 within the central region (en¬ 
closed by the dashed line), but increases rapidly outside this 
region. 

Outside the boundary, the large relative error simply fol¬ 
lows from the very small exact values |(igx|. To illustrate this 
point, we consider d{„ „{0) with \m\ - \n\ - + j. In this case, we 
have an analytical expression 


dUe)-(-^y- 


V 


1/2 


0\j+'” / 

cos-j (sin-j 


Using the symmetry, one can obtain d^^^jiff) - [- sin(0/2)]^''. 
For j - 100 and 6 = jijb, one can easily obtain the exact val¬ 
ues dy,^j(7Tl6) - 3.974 X 10“"*, which lie outside the bound¬ 
ary [see the dashed line of Fig. |2l[d)]. By contrast, although 
the numerically calculated values d^^.^.{nl6) ~ 10“" are very 
close to zero, they are much larger than the exact values, lead¬ 
ing to a large relative error Ajgi. 


(a) (b) (c) 
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As shown in Fig.[T] we plot the computed results of di,^„(0) 
in the plane (m, n), with m, n 6 {-j, + /]. For a relatively large 
spin j - 40 and a given 0, is appreciable only in the 

central region and tend to zero quickly outside this region. 
The boundary of this region is determined by (see also the 
dashed lines of Fig.|2]i 

- 2mn cos 0 = j(J +1) sin^ 0, (9) 

at which d'‘dl, „{0)/d0^ = 0 for k - 1,2. This boundary equa¬ 
tion follows from the differential equation of the d-matrix 
Similar boundary equation has been obtained using the WKB 
approximation 12^ . 


FIG. 2: (Color online) The absolute error A^bs = K'fcomp - dcx)l (the 
upper panel) and the relative error logj^ Apgi (the lower panel) of the 
Wigner’s d-matrix dm,„(6) with the spin j = 100, and the rotation 
angle 6 = n/6 (left), n/A (middle), nil (right), respectively. The 
computed results are obtained by diagonazing the matrix in Eq. 0}, 
using the ZHBEV subroutine of LAPACK; While the exact results 
are obtained by MATHEMATICA 10.0. The dashed lines are 
given by Eq. ®. 

Due to the same reason, the d-matrix elements for other 
values of 0 also show large Apei outside the central region. 
For example, 0 - nl2, the boundary of the central region 
is a circle: - j(j +1) with a radius ~ j. The d- 

matrix elements with \m\ - \n\ - j have the same exact value 
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dL,ni^l'^) - 1/2-' ~ 7.9 X 10 which is much smaller than 
the numerically calculated value ~ 10“'^, yielding a large Ajei. 
The computed „(7r/2) at mn - 0 also show large relative er¬ 
ror. To see it clearly, we use the exact result Jit] 

where PJix) is the associated Legendre polynomial. When 
0 = 7r/2, we obtain the exact results - 

(-l)”Wj^^j(;r/2) oc f’'j'(0) = 0 for odd j - m. In contrast, the 
computed results are small but nonzero, thus Arei —> o°. 



FIG. 3: (Color online) The maximum absolute error of di,^„(9) (solid 
circles) and that of ddl„ ,„(d)ld6 (open circles) as functions of spin j. 
Blue solid line: the fitting result, Aabs.max ~ + b) x 10^'“* with 

a = 2.568 X 10““' and b = 0.5758. The precision of the d-matrix can 
reach 3.275 x 10“*'* at j = 100. Red dashed line: maximum error of 
the first-order derivative ~ j x Aabs.max- 

Finally, we discuss the scaling of the error in evaluating 
the d-matrix and its various derivatives with increasing spin 
j. For this purpose, we sweep m, n, and 6 and calculate the 
maximum absolute errors Aabs.max for and dd'ln „i9)ld6 

as functions of j for j up to 100. Since the maximum absolute 
error of the d-matrix almost appears inside the boundary (see 
the upper panel of Fig. |2]i, we only sweep (m, ri) within the 
central region and increase 9 from 0 to 7r/2 with an increment 
7r/36 im. We find that typically the maximum absolute error 


appears at m ~ n and 9 ~ 0 or :7r/2. As shown in Fig. [2 one 
can find that numerical results of lO'^ x Aabs.max can be well 
fitted by af + b, with a ~ 10“"^ and b ~ 0.6 (see the solid 
line). This precision is slightly worse than the previous one 
Aabs.max ~ for the spin y upto 100 1^ . However, 

if the scaling persists to larger j, our method could provide a 
smaller error as j > 405. One can also note that the max¬ 
imum absolute error of ddi, „{9)ld9 can be approximated by 
j X Aabs.max (the dashed line). More generally, from Eqs. (IH 
and (|7]), we can deduce that the maximum absolute error in 
evaluating the ^th-order derivative is larger than 

that of the d-matrix by a factor 0{ f ). 

V. CONCLUSION 

In summary, we have presented a very simple method to 
evaluate accurately the Wigner’s d-matrix by diagonalizing 
the angular-momentum operator 7, in the eigenbasis of J^. 
The coefficients of Fourier-series expansion of the d-matrix, 
closely related to the eigenstates of Jy, are shown to be not 
larger than unity. This enable us to avoid the large-number 
problem. The absolute error of d'l„ „(0) can reach ~ 10“^"^ 
for spin j up to 100 and the relative error ~ 10 within the 
central region - 2mn cos 0 < j{ j + 1) sin^ 9, outside 

which the values of the d-matrix tends to zero quickly. As one 
of the main advantages, we show that for given j and 0, all 
the elements of Wigner’s d-matrix and their k-th order deriva¬ 
tives can be obtained by diagonalizing 7,- only once 1^ . The 
method presented here is a part of the larger class of algo¬ 
rithms based on the matrix exponents iH. As the matrix ex¬ 
ponent of Jy, the Wigner’s d-matrix could be computed even 
more accurately and efficiently since the matrix in Eq. (|6]l is 
Hermitian and tridiagonal. 
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